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A solution is given to the problem of representing the complex properties of ionospheric 
characteristics on a worldwide scale, including their diurnal variation, by numerical analysis 
of ionospheric data as they are obtained from a network of sounding stations, without prior 
hand operations. The problem is complicated by two basic difficulties: (1) the data are 
affected by noise (random fluctuation) and (2) the stations are irregularly positioned in the 
two space dimensions. The second difficulty is overcome by a general method for construct- 
ing functions orthogonal relative to the distribution of the stations. Special filtering proc- 
esses are employed for the optimum separation of noise from real physical variations. The 
end product of the analysis is a table of numerical coefficients defining a function r (X, ^, t) 
of three variables, latitude (X), longitude (6) and time (0, which can be used to compute the 
ionospheric characteristic at any desired location or instant of time. The methocl applies 
to any ionospheric characteristic; however, as a means of illustration we use in the present 
paper only the characteristic, monthly median of the /''2-layer critical frequency (foF'^). 



1. Introduction 

It is well known that the ionosphere experi- 
ences diurnal and geographic variations, as well as 
long term variations connected with seasons and the 
11 year solar cycle. Altliough many of these varia- 
tions are systematic and predictable witli a consid- 
erable degree of acciu'acy, for the /^'-layer theoretical 
models do not yet give good representation of many 
of the details. Therefore, worldwide i'epres(nitations 
of i^-layer characteristics are best macU^, using meas- 
urements obtained from a network of soiniding 
stations. 

Ionospheric maps have been made and still are 
being produced by persons having considerable 
knowledge of ionospheric data and experience in 
drawing maps. Many graphic techniques are em- 
ployed for smoothing raw data and for eliminating 
inconsistencies. However, such procedures are rela- 
tively slow and tedious and, to some extent, are not 
repeatable. The accuracy, as well as repeatability, 
are dependent upon the skill, knowledge, and expe 
rience of the people carrying out the subjective steps 
of the process. Moreover, as a means of simplifying 
hand operations, unrealistic assumptions have been 
made, such as the ' 'longitude zone system,^' in which 
no longitudinal variation is taken into account or 
represented within the defined longitude zones. 

Several years ago, the authors attacked the prob- 
lem of representing the complex properties of iono- 
spheric characteristics on a worldwide scale, including 
their diurnal variation, by numerical analysis of iono- 
spheric data as they are measured at the stations, 
without prior lumd operations.^ The problem is 



* A sequel paper on "Methods for applying numerical maps of ionospheric 
characteristics" will appear in an early issue of this journal. 

1 The need for mapping methods based on numerical methods and the use of 
high-speed computers has been felt for several years [CCIR, 1959]. 



complicated by two basic difficulties: (1) the data 
are affected by noise — random fluctuations pro(hiced 
from a number of sources (ch. 4) — and (2) the iono- 
spheric stations (fig. 1) are irr(\<>iilarly positioned on 
the earth. The noise would produce a very rough 
and physically unacceptable map if the original data 
were represented exactly without time and space 
smoothing. Thus a certain amount of smoothing is 
necessary, but too much smoothing would give a map 
which does not respect the true physical variation as 
well as possible. The irregular distrihuiioii of stations 
presents prohlems both in data-fitting (ch. 2) and in 
preserving the stability and physical soundness of 
the representation in areas where few if any stations 
are available (ch. 5). Also the fact that the set of 
available stations varies from month to month com- 
plicates the data-fitting processes and the comparison 
of numerical representations for different months. 
In the past, people experienced in drawing maps fiave 
overcome these difficulties more or less intuitively 
using empirical knowledge of the ionosphere and good 
sense. It is not easy to give these qualities to a com- 
puting machine. 

The solution^ to the problem consists of well- 
defined mathematical operations — described and 
illustrated in the present paper — which have been 
programed for use on several large-scale digital com- 
puters.^ Input to the computer program consists 
of the measmements of an ionospheric characteristic 
from all available stations for a given month. The 
diurnal variation is represented by Fourier analysis of 
the 24 hourly measurements from each available 



2 For a brief summary the reader can refer to [Jones and Gallet, I960]. 

3 The computer program has been developed completely for botli the IBM 704 
and 7090 computers and in part for the GDC 1604. 
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Figure 1. Map of ionospheric stations for December 1957 



station (sec. 3.1). Then the worldwide geographic 
variation of each Fourier coefficient is expanded in a 
series of functions analogous to surface spherical 
harmonics (sec. 3.3). The optimum separation of 
noise is obtained by truncation of orthonormal series 
(ch. 4). At the end of the analysis the diurnal and 
geographic variations of the characteristic are repre- 
sented by a relatively small table of coefficients 
defining a function of three variables: latitude, longi- 
tude, and time of day. Such a function is referred 
to as a ^ ^numerical map'' (ch. 7). 

The methods used here are general enough to be 
applied to any ionospheric characteristic. In fact, 
they have already been successfully used to represent 
such characteristics as : the critical frequency if 0^2) , 
the 3,000 km maximum usable frequency factor 
(F2-M3000), the maximum electron density (A^- 
max) , the height of A^-max, and the quarter thickness 
of a layer. With only slight modification the 
methods can also include variations with height 
above the surface of the earth. For the sake of 
illustration, we restrict ourselves here to the charac- 
teristic/0^2 monthly median, since this characteristic 
is the most important one for radio propagation, and 



its variations are the most difficult to represent. In 
chs. 2 and 6, we outline briefly the mathematical 
methods emplo3^ed in this analysis. 

2. General Data Fitting Method 

As was mentioned in the introduction, the geo- 
graphic variation of each Fourier coefficient (obtained 
from the diurnal analysis) is represented by series of 
functions analogous to surface spherical harmonics. 
The approach generally used in global analyses of 
geophysical data has been to first draw contour maps 
by hand and then to analyze in spherical harmonics 
the values read from the maps at the intersections of a 
regular grid. However, such methods incorporate 
many of the undesirable features of the hand opera- 
tions used currently to produce ionospheric maps 
(ch. 1). 

We have attacked this problem from the opposite 
direction by first analyzing the data directly as they 
are obtained from the stations. Then contour 
maps are computed /rom the analyses when desired.^ 

4 Examples of such contour maps are shown in figures 11a, lib, 12, and 13. 
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The method used is iiiathematically wcU ck^fined, 
enth'ely repeatable, well adapted to automatic com- 
puting and, in a sense, more objective than the 
hand methods previously mentioned. As a result of 
this new approach, we have had to face some addi- 
tional complications (ch. 5). In the expansion of 
the geographic variation it is natural to use 
orthogonal ^ functions, since an optimum cutoff of 
the series has to be made (ch. 4). The classical 
spherical harmonic functions [Byerly, 1893, pp. 
144-218] could not be used, however, since they are 
not orthogonal relative to the positions of stations. 
Thus we have had to construct functions — analogous 
to spherical harmonics (sec. 3.3) — which are orthogonal 
with respect to the irregularly spaced coordinates of 
the stations. The problem is even more complicated 
by the fact that the set of available stations varies 
from month to month, making it necessary to con- 
struct new orthogonal functions for the analysis of 
each month's data. 

For this purpose we have employed the very 
general data-fitting method described briefly in this 
chapter. Of course, the method also applies to the 
Fourier analysis of the equally spaced measurements 
of the diurnal variation. In the following discussion, 
the points Xi (where the measurements are made) 
can be any set 



Xf \Xi\) Xi2) 



^iMJ 



= 1,2, 



A^ 



irregularly positioned in an M dimensional space. 
The coordinate functions Gjdx), which are fitted to 
the data by the method of least squares, can be any 
linearly independent functions of M variables (e.g., 
elementary transcendental functions, Chebychev 
polynomials, or spherical harmonics). 

2.1. Least Squares Method 

Let ^1, ?/2, . . ., 7/iv denote measurements of an 
ionospheric characteristic taken at the points Xiy 
X2y . . ., Xn> Suppose that a class oj junctions 
Ya{x) is given and a criterion for evaluating how 
well each of these functions fits the values ?/,. (We 
will take Yaix) as linear sums of coordinate functions 
Gkix)) Then the problem of data-fitting is to 
determine that function Y{x) from the given class 
which ''best'' fits the yi relative to the given criterion. 
We have chosen the least squares criterion since it is 
well-adapted to the analysis of data affected by 
noise and it is far easier to compute than other 
methods such as the minimax [Stiefel, 1959]. The 
question of choosing the proper classes of functions 
Ya{x) is discussed in chapter 3. 

We assume given a set of linearly independent 
coordinate functions, Gq{x), Gi(x), . . ., Gk{x), with 
K<iN, and it is required by the least squares 
criterion to find that function 



YAr)=i:D,GM 



k=0 



(1) 



for which the sum of squares of residuals 

6 The t(>rni "ortliop;onar' is used here in the sense of a discrete distribution 
[Szcgo, 1959, pp. 33-37]. 






(2) 



has a minimum value with respect to all real coeffi- 
cients Djc.^ It is easily shown that necessary and 
sufficient conditions for Ek to be minimized are 

^=0, ^=0, 1, . . ., K. (3) 

These give the well known normal equations 

S D,(G„,G,) = iy,G^) , m=0,l,...,K (4) 

/c=0 

where we adopt the notation 



and 



N 

(GmfGjc) =S Gmi-^i) Gjc (Xf) , 
i = l 



{y,Gm)=J2yiGm{,Xi), 

i=l 



(5) 



and we will later use 



N 



iy,y)=J2yl 

i=l 

For solving the normal equations (and for a number 
of other reasons to be discussed) we have made 
extensive use of orthogonal functions. 

2.2. Gram-Schmidt Orthogonalization 

A set of functions Ao(x), Ai{x), . . ., Ak(x) is 
said to be orthogonal with respect to the points 

Xi, X2, . . ., X]sr, II 



(Aj„A^)=0 when JcT^jn. 



(6) 



Thus we see that if the coordinate functions Gjc{x) 
were orthogonal, the system of normal eq (4) would 
be uncoupled and its solution would be simply 

A=^^^. k=Q,l,...,K. (7) 

On the other hand, we can apply the Gram- 
Schmidt orthogonalization process [Davis and 
Rabinowitz, 1954]^ to form the system of functions 



Ao(:x) = Go{x) 



(8) 



Ajc {x) -2 ajcpA^ (x) +Gic (x) , ^-1,2,.. ., K, 

p=0 



6 For simplicity we have written Dk witli only one subscript k. However, we 
must remember that the terms in (1) are generally not independent and therefore 
all of the coefficients depend on the value of K. 

7 Our procedure differs from that of Davis and Rabinowitz [1954] in that we 
have reversed the order of orthogonalization and normalization in order to reduce 
the accumulative rounding error. (See sec. 2.3 and also ch. 6.) 
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satisfying (b), where the coefficients a^,, are given by 



(^kp 



{Ap,Ap) 



(9) 



Then the least squares solution Yk(x) can be written 
in the form 

Ydx)=j:a,A,(x) (10) 



where 



%= 






^-0,1, . . .,K. 



(11) 



The orthogonal series (10) is, of course, sufficient 
for many purposes, but in our applications it is 
useful to have Yk{x) in the simpler form (1). For 
this purpose we compute the triangular matrix 



^00 
^10 



^01 

^11 



h 



02 
ll2 



4 






lo 



(12) 






Ik- 



where the elements Ij,^ are obtained by 

OSk^K-p 



(13) 



starting with the top row and going from left to 
right. Then we have for the desired coefficients 



J^k—lK-k,k' 



(14) 



A number of advantages are gained from the use of 
orthogonal functions. As was previously shown the 
Aj^{x) uncouple the system of normal equations so 
that terms in the series (10) are independent. As a 
consequence we obtain automatically the least 
squares solution Y]^{x) for all degrees from zero 
through K, and therefore can test the physical 
significance of each term to determine where the 
series should be truncated (ch, 4). Moreover, the 
computation of Ek is greatly facilitated since by 
substituting (10) into (2) and applying (6) and (11), 
We obtain 

EK-^{y,y)-l:al{A,,A,), (15) 

/c = 

Hence Ej, is computed recursively for each degree by 
E,^{y,y)-al{A„A,) (16) 

E,=E,^,-al{A,,A,) for l^k^K. 

It should also be mentioned that in the construction 



of the orthogonal functions Ak{x) we use only the set 
of points Xi and the functions Gk(x). Consequently 
when several different sets of values yi are measured 
at the same set of points Xt and are to be fitted by 
sums of the same functions Gi,{x), only one orthogonal 
system is required to uncouple all of the resulting 
systems of normal equations. Thus in the repre- 
sentation of geographic variations of Fourier coeffi- 
cients, the same set of orthogonal functions can be 
used for all Fourier coefficients for a given month. 
This is an important factor in the economization of 
computer time, since — as will be shown in chapter 
4 — there are 17 Fourier coefficients to be represented 
in the analysis of /oi^2. 

2.3. Normalization 

A set of functions Fi{x), F2{x), . . ., Fk(x^ is said 
to be orthonormal with respect to the points Xi, 
X2, . . ., x^r if, in addition to the orthogonality 
conditions (6), they also satisfy 



{F,,F^) = l fori-0, 1, . . 
When Yk{x) is expressed in the form 



K. 



YK{x) = ^d,F,{x), 

k=0 



(17) 



(18) 



the orthonormal coefficients d„ have the simphfied 
formula 



and 






iy,F,) k-- 



--0,l,...,K 
l^k^K. 



(19) 
(20) 



One main advantage of this normalization is that the 
significance of each term d^Fjcix) in (18) — determined 
by the reduction it makes on the sum of squares of 
residuals (20) — is seen in the relative size of dl. 
Another advantage is obtained from the simplified 
interpretation of these functions as unit (orthogonal) 
vectors in an N dimensional vector space (ch. 6). 
Using the relation 



F,{x) = 



A,{x) 



V(^*,A) 



(21) 



and (8), the generation of the orthonormal functions 
Fjc{x) becomes 



k-l 



where 



and 



FjXx) = TjdkpF,{x)+d,,G,{x) 

p=0 



(22) 



dkp- 



--cikj 



V {Ajj,Aj,) 
{A„A,) 



OSp^k-l 



da= 



V(A,A) 



(23) 



Then from (19), (22), and (11), the coefficients in 
tfie solution (18) become 
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4=a*V(^(i,-4t). 



(24) 



The application of the Gram-Schinidt ortliog-onal- 
ization process to least squares problems has been 
treated in a paper by Davis and Rabinowitz [1954]. 
For completeness we included here an outline of 
their development with one modification^ — the order 
of ortliogonalization and normalization is reversed. 
The advantage of this modification is that accumula- 
tive rounding error can be considerably reduced. 
Whereas we compute {Afc,Ak) as a sum of squares, 
the procedure in the reference cited above is to com- 
pute it as a difference of positive numbers 

which can be so small as to result in total loss of 
significant digits. As can be seen from (23), errors 
in {Ajc,Ak) have a strong accumulative effect on the 
construction of tlu^ functions Fk(x). 

We liave also foimd, however, that a very serious 



additional rounchn 



' accumuhites in the CJ rain- 



Schmidt orthogonalization process when the coordi- 
nate functions Gfdx) are far from being orthogonal. 
That such rounding errors can have a sigiuficant 
effect on the least squares solution is shown in cluipter 
6 by a numerical example with errors occurring by 
as much as 30 percent. It is well known that similar 
types of errors occur when solving the normal cq (4) 
by other metliods.^ Thus the danger of accujnulative 
rounding error is not so much a characteristic of the 
particular method we have used, but is an inherent 
difficulty in all large-scale least squares problems. 
Also included in chapter 6 are an explanation of the 
process of accumulation of the error, a quantitative 
method for estimating its effect, and a reorthoijonaU 
izaiioth process which keeps tlie rounding error unckn- 
control no matter liow large the system. 

For the diurnal representation of th(^ (hit a, we 
use Fourier series which are particular orthonormal 
series for equally-spaced data points. For the geo- 
graphic variation we employ an orthonor-nuil series 
constructed from the procedures described above. 

3. Choice of Functions 

3.1. Diurnal Variation 

The most natural method for representing the 
diurnal variation is Fourier analysis, since iono- 
spheric characteristics are periodic functions of time. 
Moreover, the trigonometric functions associated 
with Fourier analysis are automatically orthogonal 
with respect to the equally -spaced points of measure- 
ment; hence the computational problems of least 
squares fitting are greatly simplified. Since the 
observations are made at each hour in the zone 
time (ZT) of each station, they cannot be intercom- 
pared for different stations until time corrections 

? The process of accumulation of rounding error is similar to that described by 
Lanczos [1956, pp. 123-130]. 

9 See for example [Lanczos, 1956, pp. 118-122; Kunz, 1957; Forsythe, 1957, p. ;77 
Forsythc and Rosenbloom, 1958, pp. 20-21, and references contained therein]. 



are made. As will be shown, Fourier analysis pro- 
vides a simple means for such corrections. The 
most important use of Fourier analysis, however, 
is in the separation of ''noise'' (random (MTor) from 
the ''real" diurnal variation of the data (ch. 4). 

We give here the essential formulas employed; 
for more comprehensive treatment, the reacU^* can 
refer to numerous texts on Fourier analysis, l^et 
X denote the zone time hour angle (degrees) (h^fiiuMl by 



x=15° (ZT)-180^ 



(25) 



where ZT is given in hours. Thus, for exainj)le 
x=0° at noon (ZT). We let yi, 7/2, . . ., y24 denote 
the houily measurements corresponding to the fiours 



x, = 15°i i=l,2, . . .,24, 



(26) 



respectively, and we choose 1, cos jx and sin jx 
{j=\, 2, . . ., H, and 2i?+1^24) as coordinate 
functions G'k[x). It is well known that these functions 
are orthogonal with respect to the points (2()). 
Therefore from (7) the least squares solution has 
the form 



5^2/f-M(^)=^«o+S h; (*os jx+Pj sin jx], (27) 
7=1 



where 



I 24 

o^j=j2^ y' cos^x,-, 



I 24 

^i=Ti^ Z) Vi sin jxt 



(28) 



^^j^ll- 



Equation (27) can also be written in tlie convenient 
form 



J^2H+i(2)=ao+S Cj t'OS iJJ--<t>j) 



(29) 



4'i=<t>j+j(e-e^), 



(31) 



6 is the longitude (degrees east of Greenwich) of the 
station. Therefore the representation of the diurnal 
variation takes the form 



where the amplitude Cj and j)hase (jyj are given by 
Cj=^|a)+^ and </)y=arctan -- (30) 

Corrections to Local Mean Time 

In order to intercompare data from differc^nt | 
stations it is necessary to correct for the small 
difference between the actual local mean time (LMT) 
at each station and the time at the reference longi- 
tude dji of the zone. Sucli corrections can be made 
very simply in (29) by a shift of the ])li{ise wJiere 
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// 
U2H+i{t) =ao+X) [(ij COS jt + bj sin jt] (32) 
;=i 

where the Fourier coefficients corrected to LMT are 
given by 

aj--Cj cos \l/j bj=Cj sin \l/jy (33) 

and t denotes the local mean hour angle, 

^=15°(LMT)-180°. 

3.2. Main Latitudinal Trend 

The Fourier coefficients are corrected to LMT 
(sec. 3.1) so that their main geographic variation 
becomes latitudinal and is therefore greatly simpli- 
fied.^^ Among the first problems to be solved for 
representing this variation was that of selecting a 
suitable set of coordinate functions Gjc(x) for the least 
squares fit. Polynomials seemed to be the most nat- 
ural type of function, but the question arose as to 
what would be the best independent variable for the 
polynomials. The simplest variable tried was the 
geographic latitude X. However, when a sufficiently 
high degree was taken, the polynomials in X became 
unstable (i.e., wildly fluctuating) in regions such as 
near the poles where little or no data were available 
(figs. 2 and 3). Much more stable representations 
were obtained by using polynomials in sin X. More- 
over, it was found that the equatorial variation 
could be represented in more detail by these functions 
than by the polynomials in X of the same degrees 
(figs. 2 and 3). 

10 An atlas of graphs of the geographic variations of the Fourier coefficients for 
/0F2 median for four seasonal months has been prepared and will soon appear as 
an NBS Technical Note [Jones, 1962]. This atlas illustrates the very system- 
atic aid well-defined variations of these coefficients. 
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Figure 2. Representation of main latitudinal variation of 
Fourier {tirne series) coefficients slq for U'P-i monthly median 
by least squares polynomials of degree 10. 
December 1957 96 stations. 



One explanation of the superiority of the variable 
sin X is the following. The fitting of polynomials in 
sin X to data located at a set of latitudes Xj is equiva- 
lent to fitting polynomials in a variable x to the same 
data located at the correspondingly shifted set of 
points Xf=sin X^-. The shifting of the data resulting 
from this transformation has the effect of pulling the 
data symmetrically away from the equator toward 
the poles. Thus the data become more uniformly 
distributed in the interval — l^a^^l, and the 
sharpness of the variation near the equator is reduced. 
To illustrate this spreading effect the same data 
shown in figures 2 and 3 have been plotted against 
sin X in figure 4, together with the polynomial of 
degree 10 in x=sin X. 

A second explanation of the improved behavior of 
the polynomials in x=sin X can be given in terms of 
the corresponding orthonormal functions. Choosing 
powers of the independent variable as coordinate 
functions, we generate (sees. 2.3 and 2.4) one ortho- 
normal system Fj^^^ (X) corresponding to the station 
latitudes X^ and a different system Fi^^ (sin X) cor- 
responding to the shifted set of points a::i=sin X^. 
The resulting least squares representations take the 
forms 



n'o'(x)=i]4iiFi»(x), 



(34) 



Yi'^ {sin X)=T. dp F^' 



(2) 



(sin X), 



respectivel}^. A comparison of the different ortho- 
normal functions (for degrees 8 to 10) is shown in 
figure 5. To simplify the comparison, an additional 
normalization was made so that all of the graphs 
have a common value at the right end point. It can 




Figure 3. Representation of inain latitudinal variation of 
Fourier {time series) coefficients ao for foF2 monthly median 
by least squares polynomials of degree 13, 
December 1957 96 stations. 
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Figure 4. Representation of main latitudinal variation of 
Fourier (time series) coefficients ao for foF2 monthly median 
by least squares polynomial {degree 10) in sin X. 
December 1057 06 stations. 



be seen that at high degree (8^^^ 10) the poly- 
nomials in X have very large values near the right 
end point compared with their values near the 
center. Thus to represent the strong geographic 
variation near the equator the polynomials F\^^ (X) 
of high degree have to be multiplied by large coef- 
ficients d]^\ which results in the blowing up effect 
near the poles (figs. 2 and 3). In contrast, for the 
polynomials F^^ (sin X) the maximum ainphtude 
of oscillation is more uniform throughout the int erval, 
hence the greater stability of Y)^^ (sin X) near the 
poles. In addition we see that near the center of 
the interval the distances between successive maxima 
and minima are less for the Ff ^ (sin X) than for the 
Fj^^ (X). Thus we can understand the better rep- 
resentation of equatorial variation by Yf^ (sin X). 

To conclude this discussion we note that the 
behavior of the polynomials F[P (X) resembles that 
of the classical Legendre polynomials [Byerly, 1893, 
pp. 184-185]. On the other hand, the amplitude 
distribution of the F^^^ (sin X) more nearly approx- 
imates the exactly uniform distribution of the 
classical Chehychev polynomials [Jones, Miller, Conn, 
and Pankhurst, 1946, pp. 194-195]. It is well known 
that orthogonal series of Legendre ])olynomials tend 
to blow up TU^ar the end points [Lanczos, 1938, pp. 
144-145], whereas series of Chebychev polynomials 
minimize the maximum vvvov, and therefore are as 
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stable at the end points as they are at the center. 
Thus by analogy one would expect the behavior 
demonstrated in figures 2 and 3. 

A similar study was also made witli polynomials in 
the variable sin [90° sin X], which is the natural 
extension of the spreadhig process described above. 
In this case, the representation was improved even 
more around the equator, but significant geographic 
variation in the temperate latitudes was squeezed 
into the poles and so lost. Therefore polynomials in 
sin X were chosen for representing the main latitudi- 
nal trend. The determination of the ^'best degree" 
for these polynomials is discussed in chapters 4 and 5. 

3.3. Mixed Latitudinal and Longitudinal Variations 

The existence of systematic longitudinal varia- 
tion ^1 is illustrated by the graph in figure 6 of Fourier 
(time series) coefficients ai — first harmonic, cosine 
part — for /0F2 monthly median, plotted against 
latitude X. Also shown in the figure is the rep- 
resentation of the main latitudinal trend by 
means of a polynomial of degree 10 in sin X. By 
means of special plotting symbols used to signify 
approximate station longitudes (see legend), we see 
that the coefficients are not randomly dispersed 
about the main latitudirud trend, but — in certain 
regions, particularly around the equator — are sys- 
tematically arranged according to longitude 6. We 

11 That the longitudinal variation, illustrated in figure 6, is consistent for differ- 
ent seasons and for periods of high and low solar activity is demonstrated by the 
atlas of graphs of Fourier coefficients for /0F2 monthly median referred to at the 
beginning of section 3.2. 
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discuss here the choice of functions for representing 
these second order mixed latitudinal and longitudinal 
variations. 

The usual method of representing global variations 
of geophysical phenomena is by spherical harmonic 
analysis [Byerly, 1893, pp. 144-218]. However, the 
classical surface spherical harmonics — satisfying La- 
Place's equation — are not orthogonal with respect 
to the positions {\i,di) of the ionospheric stations. 
Therefore for convenience, and with no loss of 
generality, we choose as our coordinate functions the 
simplest set G/,(\,d) (see table 1) of which the surface 
spherical harmonics are linear combinations. A 
particular set of the functions G,c (^^^n be specified 
by assigning the values go, (lu and qo which are, re- 
spectively, the highest powers of sin X for terms 
involving: (1) no longitudinal variation, (2) first 
order longitudinal variation, and (3) second order 
longitudinal variation. Equivalently, we could also 
specify the values ko, ki, and k2=K taken by the 
index k at the end of each of the three groups of the 
Gi,{\fi). These values are related to the g-values 
by the equations 

h--q.o kr = kr-i + 2(qr + l) for r=l, 2. (35) 

From table 1 we see that the first group of functions 
(powers of sin X) is the same set chosen (sec. 3.2) 
for the main latitudinal trend. We note that the 
zonal harmonics are linear combinations of these 
functions, and the first and second order sectorial 
and tesseral harmonics are linear combinations of 
the first and second order terms in longitude (see 
table 1) used for the mixed latitudinal and longi- 
tudinal variation. As a consequence, a least sc^uares 
representation in the form of a linear combination 
of the G}c{^,d) has many properties of a spherical 
harmonic analysis. For example, it is periodic in 
longitude and constant at the poles, and the first 
and second order terms in longitude are weighted 
according to latitude by the functions cos X and 
cos^ X, respectively. In fact, a series in the Gi{\M) 
is identical — but expressed in a different form — to 
that which would be obtained from spherical har- 
monic analysis. 





Table 1. 


Geographic fum 


itions Gi 


(\,d) 


Main latitudinal 


Mixed latitudinal and longitudinal variation 


variation 


First order in longitude 


Second order in longitude 


k 


ak{\,d) 


k 


r^k{\e) 


k 


rrk{\d) 



1 
2 


1 

sin X 
sin2 X 


ko-\-2 
A'o+3 
fco+4 


cos X cos 9 

COS X sin d 

sin X cos X cos 

sin X cos X sin d 


fci-fl 
fci+2 
fci-f-3 
/ci-h4 


cos2 X cos 20 

cos2 X sin 20 

sin X eos2 X cos 20 

sin X cos2 X sin 20 








h 


sin ^0 X 












ki-1 
kx 


sin'i X cos X cos 
sin^i X cos X sin d 


K-l 
K 


sin52 X cos2 X cos 20 
81^2 X cos2 X sin 20 



The geographic variation of a Fourier (time series) 
coefficient is therefore represented by a function of 
the form 
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r^-(x,e)= V, dj'\{\fi)^f:. D,GA\e) (36) 

where the Ffc.(\,6) are oi'thonormal funclions (Un- 
filled with respect to the positions of the iouospherie 
stations (taking as coordinate functions the Gjc(\,d) 
in table 1), and the coefficients dt and Djc are ob- 
tained by the method of least squares (ch. 2). The 
determination of the ''best'' representation with these 
functions is discussed in chapters 4 and 5. 

4. Optimum Separation of Noise From Real 
Physical Variation 

It was previously mentioned that ionospheric 
data are affected by noise (random fluctuation) 
produced from a number of sources. The noise 
is due in part to limitations of equipment at various 
stations and to errors of scaling and rounding. 
A large part of the noise is the result of statistical 
fluctuations in the sample medians, the sample size 
being at most 31 and frequently much less. These 
fluctuations are, to some degree, caused l)y intrinsic 
variation of the physical plienoniena being nu^asured. 

The noise is evidenced by such occurrences as 
unusual roughness or unrealistic flattening of the 
diurnal plots during certain hours. The presence 
of noise is also suggested by random inconsistencies 
of data at groups of neighboring stations as observed 
from data-comparisons on a worldwide basis. In 
the present chapter, however, we give more objective 
evidence of the noise by means of statistical and 
mathematical data-analyses. Included are methods 
for quantitative estimation of the noise and for its 
optimum separation from real physical variation. 

We sliall therefore consider ionospheric data as 
sampled values yt from a function/(j') which is the sum 
of two components: (1) the real physical character- 
istic Ji(x), and (2) a random noise component 
'f2(x). Thus we consider each value y,- as a sum 



where 



yi=yi 



-Vi 



yl'^=Mx,) ^nd y?^=Ux,) 



(37) 
(38) 



The values of the noise ?/f^ are assumed to be, 
independently, normally distributed with mean 
zero and standard deviation a, a quantity taken 
to be the measure of the noise. 

Generally the noise is small compared to the main 
ph\^sical variation, but its effect must be carefully 
studied since we wish to represent the physical 
characteristic with as much detail as possible. 
The noise would produce a very rough and physically 
unacceptable representation if the original data 
were fitted exactly without time and space smooth- 
ing. Thus a certain amount of smoothing is neces- 
sary, but too much smoothing would produce a 
representation which does not respect well enough 
the true physical variation of the data. Our prob- 
lem is therefore to represent as accurately as 
possible the component /i (a:), given only the sampled 
values yt ol f(x). We treat first the separation of 
noise from the diurnal variation. 



4.1. Separation of Noise From the Diurnal Variation 

We employ a type of jnathemalical iiller ^^ wliich 
rejects that part of the ^'signal" (ionosplieric data) 
produced mainly by noise /2(;r) and accepts tiie part 
representing mostly true physical variation fi(x). 
Our main tool is the Fourier analysis of the sampled 
values yi (sec. 3.1) which decomposes the diurnal 
variation into eleven harmonics 



Cj cos {jx—\l/j)=aj cos jx-\-bj sinjx 



(39) 



and thereby gives a discrete power spectrum Cj. It is 
shown that some of these harmonics represent mostly 
j'i(x), whereas the others are produced mainly by/2(i^). 
The proper separation of these harmonics and 
resulting truncation of the Fourier series give the 
optimum smoothing (or filtering) as well as the 
desired diurnal representation. For determining 
the proper separation of harmonics ^^ we make use of 
certain properties of the Fourier spectrum which 
characterize the two component functions of f(x). 
As a by-product we also obtain a quantitative esti- 
nuite of the noise a, which can be compared with 
results from an inch^pcMuhMit nu^thod. 

a. Spectrum for a Real Ionospheric Characteristic /i(x) 

It is well known tliat a smooth, continuous function 
fi(x) of period 27r can be expanded in a Fourier 
series — of functions 1, cos jx and sin jx (j=lj 2, 
3, . . . ) — and that the coefficients, given by integral 
formulas analogous to (2(S), approach zero as j 
increases, at a rate depending upon the smoothness 
of J}{x). It can be shown, ^^ for example, that if 
fi'\x) is piecewise continuous in (— 7r,x) then the 
coefficients approach zero at least as fast as j~^, 
whereas U J^^'\x) is piecewise continuous the rate is 
j~'\ In general the size of a in the ;/"" law increases 
or decreases with the smoothness of /i(x), so that in 
a sense the smoothness of a function is characterized 
by the value of a. 

These same laws apply (approximately) to the 
coefficients aj and bj obtained from Fomier analysis 
of a sample, yi^\ yi^\ . . . , y2l\ of iiourly values 
of fi(x). Thus the squared amplitude c5 = a5 + 6? 
decreases at least as fast as j'^"^. When c? is plotted 
against j in a log-log scale, the curve defined by 
smoothing the points on the graph will have a slope 
less than or equal to —2a. We refer to this value 
as the slope oj the Fourier spectrum. It is safe to 
assume that the real diurnal variation of an iono- 
spheric characteristic has at least a piecewise con- 
tinuous first derivative. Thus the slope of its spec- 
trum should be —2 or less. 

b. Spectrum for a Random Noise Component /2 (x) 

Let^P\ y^^\ . . ., ?/2? denote a sample of values of 
the random noise component /2(^) — each value hav- 

12 This type of filter is analogous to a "low-pass electrical filter— i.e., a filter 
designed to pass only low frequencies while eliminating all frequencies above a 
certain point [Holloway, 1958]. Such a filtering process is frequently referred to 
as "smoothing." 

13 The central idea employed is suggested by Lanczos [1956, pp. 331-344]. 
H See [Jackson, 1957, pp. 1-22]. 



427 



ing an independent normal distribution with mean 
zero and variance o-^ — and consider a Fourier analysis 
of these values. The statistical distributions of the 
Fourier components (a^-, bj, Cj, and \f/j) are obtained 
as a result of the asymptotic solution for the random 
walk [Chapman and Bartels, 1940, pp. 572-582]. 
Thus it follows from a theorem of Markoff that the 
coefficients aj and bj have independent normal dis- 
tributions with mean zero and variance 0-7 12. From 
this it is shown that the phase \l/j is uniformly dis- 
tributed in —180°^rpj^ 180° and that the amplitude 

Cj has a Rayleigh distribution with mean (rV7r/24, so 
that the squared amphtude c^- has the expected value 



EicJ)' 



s — « 

6 



(40) 



Since these results are independent of the harmonic 
j, it follows that c^ oscillates about a constant value 
0-^6, the noise level, and that the Fourier spectrum 
for random noise has a slope of zero. 

c. Spectrum for Ionospheric Data 

Since the noise is small compared to the main 
amplitude of the diurnal variation, we would expect — 
from the theory outlined above — that the spectrum 
for actual data would decay rapidly for the first 
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hariiioiiics and would tlion level off around a constant 
value (7-/6. The point at which the leveling* off takes 
place separates the harmonics representing mostly 
true physical variation from those produced mainly 
by noise. Thus we determine the optimum cutofr for 
the Fourier series. It is inevitable that a few har- 
monics will be in a "twilight zone" where the ampli- 
tude of the real physical variation is the same as that 
of the noise. However, it is not of vital hnportance 
where the cutoff is made in this region, except that 
the series should not be extended too far. 

To illustrate the method we consider the monthly 
medians of /o^2 from three typical stations (fig. 7). 
These examples illustrate the relation between 
smoothness of data and slope of the spectrum. For 
Victoria, with the smoothest data, the slope of the 
spectrum (for the lower order harmonics) appears to 
be —4; for Yellowknife it is aroimd —3; and for 
Dikson Island it is approximately —2, in beautiful 
agreement with the theory. In each case the spectra 
show a tendency to level off around harmonic 7, 
indicating the effect of noise on the higher harmonics 
and the position of the twilight zone for the optinunn 
cutoff. 

The theory has also been applied to determine the 
average oj^tirnum cutojf^ using the nuMin spectrum for 
all stations for the given month (fig. 8). As was 
expected the mean spectrum is nuich smoother than 
the spectra for the individual stations, so that the 
slope and optimum cutoff are more clearly defined. 
There is a definite change in the character of the 
spectrum at harmonic 8. It is clear that harmonic 7 
is above the noise level and so should be retained, 
but harmonic S — being in tlie twilight zone — could 
either be retained or not. We have terminated the 
series after the 8th harmonic. The eft'ect of tliis 
smoothing process is illustrated for the tliree stations 
in figure 7 by the solid lines representing the diiunal 
variation. It has been found that for other months 
and other characteristics the mean spectrum is 
generally as smooth as the one shown in figure 8, 
but the twilight zone is sometimes more extended. 

A number of other stuches have been nuide to 
determine the average optimum cutoff, in addition 
to the work on Fourier spectra. From the theory 
of analysis of random noise, the phase yj/j for noise 
harmonics is uniformly distributed; thus all values 
of the phase are equally likely to occur. By com- 
paring the phase angles from all stations for a given 
harmonic (by means of polar plots of amplitude 
and phase), we have found systematic variations 
with geographic position for the lower harmonics 
and apparently random (uniform) distributions 
for higher harmonics, in good agreement with 
the previous results. Similar investigations have 
been made to determine the distributions of the 
other components, ay, 6;, and Cj. Although reason- 
ably good agreement with the previous studies was 
found, the results for these cases w^ere not so well 
defined. Tliis was to be expected, however, since 
the parameters defining the normal and Rayleigh 
distributions of tliese componcMits ar(^ suhject to 
geographic variation. 
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d. Estimation of the Noise a- 

The theory also provides a means for com])uting 
the noise <r. Taking the mean squared ami)litnde 
for the noise harmonics 

6 j=g 

as an estimate of E{c^), we use (40) to obtain a. 
Values of the noise are given in table 2 for tlie three 
sets of data illustrated in figure 7. In a similar 
manner we compute the average noise in foF2 medians 
for December 1957 



^<IPT 



=0.25 Mc/s, 



(41) 



where c^ denotes the mean squared amplitude of 
the jth harmonic from 113 stations. Calculations 
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of the average noise in /0F2 medians made over a 
period of several years have been found to be closely 
correlated with solar activity. The details of these 
studies, however, will be given in a subsequent 
paper. 

Table 2. Noise in foF2 medians for December 1957 



Station 


Noise a 


Victoria, B.C . 


Mc/s 
0.18 


Yellowknife 

Dikson Island ._ 


.28 
.58 







The same estimates of noise have been made by 
an entirely independent method. Making use of 
the distributions of daily measurements, we have 
computed the standard deviations of the sample 
medians. The good agreement of these independent 
results gives added strength not only to the esti- 
mates of noise, but also to the determination of the 
optimum smoothing. 

Each of the Fourier coefficients is equally affected 
by noise, but we have shown that for the lower 
harmonics (j^8) the physical variation is the 
dominant part, whereas for higher harmonics mostly 
noise is present. By truncating the Fourier series 
we have eliminated about 30 percent of the noise — 
the part contained in the truncated terms — at no 
loss of real physical variation. A large part of the 
remaining 70 percent of the noise is filtered in a 
similar manner by analysis of the geographic varia- 
tion of the Fourier coefficients. This problem is 
treated in the following section. 

4.2. Separation of Noise From Geographic 
Variation 

The noise is separated from the real geographic 
variation of Fourier coefficients by a ffltering process ^^ 
similar to that used in the preceding section. Having 
expanded the geographic variation in a series (36) of 
orthonormal functions FjcO^,d) as in section 3.3, we 
obtain the smoothing (or filtering) by truncating the 
series. In this case, however, we do not have the 
elegant theory associated with Fourier analysis for 
determining the optimum cutoff. 

The metliod employed is based on the residuals 
between the original Fourier coefficients being fitted 
and the corresponding values of the representation 
function F/^(X,^) given by (36). Since orthonormal 
functions were used in the development of Yk(\,6), 
we can inspect the residuals remaining after each 
term djcFjc{\,d) is added to the series. The residuals 
approach zero as the number i^+1 of terms in the 
series is increased, and they would actually attain 
this limit (except for rounding error) when K-{- 1 is 
equal to the number of stations. However, we 
know that the Fourier coefficients are affected by 
noise; hence zero residuals are not desired. 

The criterion adopted for determining the optimum 
cutoff is to minimize the standard deviation of the 

15 See the discussion on "smoothing of observational data by the method of 
least squares" [U.S. Dept. of Commerce, NBS, AMS 9, 1952, pp. 16-18]. 



residuals. We take as an unbiased estimate oj the 
variance oj the residuals 



el- 



E, 



N-k-l 



(42) 



where Et^ is the sum of squares of residuals (20), A^ 
is the number of stations, and {N—k—\) is the 
number of degrees of freedom remaining after sub- 
traction of one degree for each term in the series 
Y^{\e) [Kendall, 1951, pp. 59-61]. From (42) it is 
possible for the estimate of the variance el (or of the 
standard deviation e^) of the residuals to increase as 
k increases, since the diminishing of E^ may become 
very slow after a certain degree. For simplicity we 
shall refer to e^ as the standard deviation oj the 
residuals. In the analyses under present considera- 
tion, the values e^, usually decrease quite rapidly at 
first, then taper off and cease to diminish appreciably 
after a certain point. This point defines the desired 
optimum cutoff. 

In order to determine this cutoff objectively we 
make use of a statistical test for the significance of 
each coefficient d^, (in (36)) based on the ''Student's'^ 
t distribution [Fisher, 1922]. Under the usual 
assumption that the observed values yi have inde- 
pendent normal distributions about a mean regression 
surface with common variance o-^, it follows that 
the quantity 



where hj^ is the expected value of d^t, has a '^Student^s'' 
t distribution with {N—k—\) degrees of freedom. 
Therefore we test the null hypothesis H^ (i.e., b],= 0) 
using a 5 percent rejection criterion 



/^^O.O 



In most of our applications the value of (iV— k— 1) 
has been large enough so that ^o.os is approximately 2. 

Although the test described above is probably the 
best possible for the present problem, we note here 
the following limitation. It was assumed that the 
distributions of all the observed values yi have a 
common variance o-^. It is known, however, that o-^ 
varies significantly with geographic location, so that 
the theory applies only approximately to our problem. 
As a result the representations are somewhat over- 
smoothed in certain regions and undersmoothed in 
others. Since we prefer a slightly undersmoothed 
representation (in order to represent the physical 
details as well as possible) we have chosen a 5 percent 
rejection as opposed to a 1 percent. 

As an illustration we consider the determination 
of the optimum cutoff of the main latitudinal terms 
(table 1) for representing the Fourier coefficient 
shown in figure 6. The last term found to be signifi- 
cant was 6?io and we have therefore terminated the 
series at this point. A graph of e^ through degree 15 
is shown in figure 9 and some of the values e^ and d^ 
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are given in table 3. In a similar manner we deter- 
mine the optimum ciitofFs for the first and second 
order terms in longitude (table 1). In figure 10 is 
shown the graph of ej, corresponding to the set of 
coordinate functions Gk{\, d) specified by ([0=10, 
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(li = 12, and 3^2 = 6 (sec. 3.3). From our test we find 
the optimum cutoffs specified by (i\ = 7 and ([0=5. 
In addition to statistical tests we have also been 
guided by certain physical properties known to exist 
in the ionosphere and by some knowledge of the 
morphology of our mathematical functions. For ex- 
ample, for a polynomial of degree k in sin X the (Us- 
tance between two adjacent maxhna will have to be 
at least (360/A:) degrees. Therefore, in order to 
represent the sharp dip in critical frequencies known 
to exist near the magnetic equator (figs. 2 and 3), 
we know a priori that very high (k^gree polynomials 
are needed. 

Table 3. Least squares fitting of spherical harmonic functions 
Gk(X,^) to Fourier {time series) coefficients ai for foF2 monthly 
medians 

l)(H'enib('r 1957 113 stations. 
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0.754A4277E 02 


2 C 


or 




6 


0.73723422E 


00 


-0.22231647E 01 


0.24343033E 03 


< 


< 
> 




7 


0.73989519E 


00 


-0.36166392E-00 


-0.98828048E 02 






6 


0.7433e301E 


00 


-0.97018918E-01 


-0.26629260E 03 








9 


0.7'»A37083E 


00 


0.633A3991E 00 


0.39049190E 02 








10 


0.71160986E 


00 


0.23279645E 01 


0.10270473E 03 








11 


0.71333002E 


00 


0.60874979E 00 


0.80027716E 00 








12 


0.71428508E 


00 


-0.61032332E 00 


-0.68859721E 00 








13 


0.7178829AE 


00 


0.89986758E-02 


0.34758507E-00 




q: 


LiJ 


lA 


0.71730649E 


00 


0.77229120E 00 


-0.15318519E-00 


Q Z 


ir 


O 


15 


0.72056954E 


00 


0.2A375897E-00 


-0.11218609E 02 


, < 




16 


0.57861927E 


00 


O.A2689060E 01 


-0.42479381E 01 


o 


(0 


17 


0.58156729E 


00 


0.99377523E-01 


-0.42751712E 01 


^ ^ 


y_ 


2 
O 


18 


0.5829'!»996E 


00 


-0.43228252E-00 


0.98574758E-01 


2 ^ 


CO 


19 


0.57096788E 


00 


-0.12750213E 01 


0.30209412E 02 


o ^ 


u. 




20 


0.57315171E 


00 


0.31006010E-00 


0.21029317E 02 


H i 

P 2 




21 


0.566325A3E 


00 


0.10180058E 01 


0.53866997E 01 






22 


0.569142A3E 


00 


-0.18124516E-00 


0.29883151E-00 


< Q 






23 


0.56555211E 


00 


-0,82855405E 00 


-0.22921705E 02 


s 1 






24 


0.544276A5E 


00 


-0.154B462AE 01 


-0.16010904E 02 


cr 


n 


25 


0.54640091E 


00 


0.3076A920E-00 


-0.87875202E 00 


? M 


g 


2r 


26 


0.51088715E 


00 


0.18782187E 01 


0.53212A79E 00 




27 


0.48734745E- 


-00 


0.15027802E 01 


0.90155041E 00 




^ 


28 


0,48894281E- 


-00 


-0,32661115E-00 


-0.75395824E 00 




o 


29 


0.46874104E- 


-00 


0.13582557E 01 


0.18835551E 01 




? 


30 


0.46133326E- 


-00 


0.88582236E 00 


0.19063634E 01 



5. Stability of the Geographic 
Representation 

To obtain a ''best" representation of the geographic 
variation, one must consider more than the optimum 
separation of noise (filtering) — i.e., more than just 
the residuals at the stations where data are given. 
The heavy grouping of stations in some regions such 
as Europe (fig. 1) and the absence of stations in other 
regions, particularly in the oceans and near the poles, 
tends to produce a sort of ''mathematical instability'' 
in the representation function — that is, unrealistic 
behavior in the areas where no data are available. 
This behavior is somewhat analogous to the large 
fl.uctuations which arise in Lagrange (polynomial) 
interpolation of clustered values with small vaiiation 
[Lanczos, 1956, pp. 346-358]. 

As we pointed out in chapter 2, the usual approach 
in the global analysis of geophysical data has been to 
first draw contour maps by hand (using the actual 
data, empirical knowledge and experience) and then 
to analyze in spherical harmonics the values read 
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from the maps at the intersections of a regular grid. 
Thus the (synthetic) values used in the analysis are 
uniformly spaced and the mesh of the grid can be 
made sufficiently small to prevent instabilities in the 
representation function. This method has the dis- 
advantage of introducing error (noise) both in draw- 
ing and in reading the maps. Moreover, the hand 
work involved is relatively slow and cumbersome. 
We have eliminated these two difficulties by first 
analyzing the data directly as they are obtained at 
the ionospheric stations. However, as a result we 
have had to face still another problem, that of 
mathematical instability described above. 

For a region where there are no available stations, 
the representation function gives a two-dimensional 
interpolation determined by the data from all 
stations, both near the region and far away. Thus 
the best representation that could be expected is a 
smooth continuation of the variations from surround- 
ing stations. Strong departures from such smooth 
continuations are probably artificial and physically 
meaningless. 

Examples ^^ of instability can be seen in the contour 
map (fig. 11a) computed from an analysis using the 
coordinate functions Gtc{\ 0) (table 1) specified by 
qo=10, qi = 7, ^2=5 — these were the optimum 
cutoffs (sec. 4.2) from the viewpoint of separation of 
noise. Included in the map are the residuals 
between the original Fourier coefficients and the 
corresponding values of Y^sCX, 6). Although the 
fluctuations in this map are not large, there is a 
definite appearance of instability in the southern 
hemisphere J particularly near the pole. By further 
truncating the orthonormal series Y^sO^, 6) so as to 
correspond to the set of ^^^(X, 6) specified by (Ig=10, 
^1 = 6 and g2=2, very little damage is caused in terms 
of increased residuals. However, as can be verified 
from the map in figure lib, considerable advantage 
is gained in terms of increased stability near the 
south pole. Table 3 gives the values of ejc and the 
coefficients djc and Djc corresponding to these functions. 

Thus it is sometimes necessary to make the cutoff 
slightly above the optimum level (determined by 
noise) in order to preserve the stability and physical 
soundness of the representation function. The 
amount of increase required in the residuals depends 
upon the world wide coverage of ionospheric stations. 



6. Reorthogonalization and the Accumula- 
tion of Error 

6.1. Accumulation of Error in the Gram-Schmidt 
Orthogonalization Process 

As was pointed out in chapter 2, the normal 
equations 



S a,{AkAm) = {y,Am) m=0, 1, . . ,,K (43) 

k=0 



(analogous to (4)) are uncoupled by the condition of 
orthogonality (6) imposed on the Ajc(x). Thus all 
terms on the left side of (43) vanish except those of 
the form am{Am,Am), so that the coefficients ajc are 
given by (11). Although (11) is algebraically exact, 
it will be shown that the coefiicients thus obtained 
could be strongly affected by accumulative errors of 
rounding. We let al^ denote the theoretically exact 
solution to (43) and write for the error in ajc (as 
computed from (11)) 



daf. — ajc ajc Zu <^mAmk 

771 = 



{Am^Aj^ 



where 

As an approximation to ^a^ we compute 

K 

Aafc= Z^ (J^mArak' 
OT = 

rriT^k 



(44) 
(45) 

(46) 



Thus if some of the niunbers AmTc, ni^k, differ sig- 
nificantly from zero (i.e., have nonzero digits among 
those retained by the computer) , it is likely that large 
values of Aa^ will occur. When we look closely at 
the orthogonalization process (sec. 2.2), it is not too 
surprising to find that some of the Amk are signifi- 
cantly different from zero, particularly if the coordi- 
nate functions Gjc{x) are not approximately ortho- 
gonal. This can be seen as follows. 

The functions Ajc{x) are constructed successively 
according to (8) and (9) . In fact, multiplying (8) by 
Am(x), m<ik, and summing over the Xi gives 



k-l 



(Aj„AJ==J2 a,^(A^,Am) + (G,,Am)^ (47) 

Then imposing the orthogonality condition (6) on the 
left member above, we have 



P9^m 



(48) 



and (9) follows by neglecting the terms —ajcpApmy 
{pT^m) which are (theoretically) zero. From (48) 
we obtain as an estimate of the relative error in a^m 



A/7 ^-1 A 

— -^ • ^kp 
^km p = (^km 



-IZiGk^Fy^f^ (49) 

p=0 \^k,J^m) 

P9^m 



16 Other examples in one dimension are illustrated in the behavior near the 
poles of the polynomials in X (figs. 2 and 3.) 



making use of (9), (21) , (23) , and (45). At this point 
it is convenient to consider the functions G]c(x), 
A]c{x), and Fjc(x) as N dimensional vectors whose 
components are the functional values at the points 
Xi, X2, . . ., x^. Then a number such as {G}c,Fp) 
can be interpreted as the component of the vector 
Gic in the direction of Fp. For the present we shall 
consider a typical term on the far right side of (49) ; 
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^x LONGITUDE 
150° 180" 210° 



240" 270* 




FUNCTIONS q = 10. q, = 7, q -- 5 




FUNCTIONS: q^'IO, q, -S.qg'Z 

Figure 11. Maps of Fourier (time series) coefficient ai (Mc/s) and residuals from analyses made with geographic functions 

Gk (X, 6) for foF2 medians. 

December 1957 113 stations. 
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thus we can assume, without loss of generahty, that 
the Gu are normahzecl in the sense of (17). Suppose 
now that G^^, has a large component in the direction of 
Fj, (say (Gk,Fj,) = %) but a very small component in 
the direction of F„„ so that {G^^Fm) is of the same 
order of magnitude as (Fp,Fp). Then the resulting 
error in a„rn could be as large as 50 percent. 

It is clear from the preceding discussion how 
significant errors can be introduced, the most import- 
ant factor being the direction of each Gf,-. relative to 
the preceding Fj„ p<ik. Although such an error 
may be small, its effect accumulates rapidly since 
each subsequent Ajc{x) is constructed in terms of all 
preceding Am{x), m<ik. Thus a small initial error 
in one Au{x) can have disastrous effects on the later 
steps of the process. 

As an illustration we consider the geographic 
representation of the Fourier coefficient, used in 
section 4.2, expanded in terms of the coordinate 
functions G],(\,6) specified by go=10, (Zi^^G, ^2=2 
(see table 1). The least squares coefficients a^c 
computed from (11) and the estimates of the relative 
error (Aajak) from (46) are given in table 4. As 
can be seen the precision of the aic for k^7 is very 
doubtful. It is shown in the following section that 
as is off by 5 percent, and other coefficients by 
considerably more, thus illustrating the usefulness 
of error estimates. 



6.2. Gram-Schmidt Reorthogonalization 

In section 6.1 we arrived at the following three 
important results: (1) large errors can be produced 
by the process of Gram-Schmidt orthogonalization, 
(2) eq (46) gives an eff'ective method for estimating 
such errors, and (3) once started, an error has a 
strong accumulative effect, but its main source 
initially is closely related to the "directions' ' of the 
coordinate functions Gk(x). If the G^^x) are nearly 
linearly dependent, the error would develop quickly, 
whereas if they are approximately orthogonal, the 
source of error would be greatly reduced. Upon 
this fact we justify a Gram-Schmidt reorthogonaliza- 
tion ^^ process, following the same lines as section 
2.2, but taking as coordinate functions the Ajc(x) 
instead of the Gj,{x). For, although the Ak(x) may 
not be sufficiently orthogonal for our purpose, they 
would in general be considerably more so than the 
original coordinate functions. Thus the initial 
source of the error would be reduced. The reortho- 
gonalization process can, of course, be repeated as 
many times as necessary to keep the error under 
control. Before giving a numerical example, we 
outline briefly the steps involved in the process. 

Following the procedure of section 2.2, we form a 
new set of orthogonal functions 

17 The notion of reorthogonalization was suggested by the Numerical Analysis 
Section (11.01) of the National Bureau of Standards. 



Table 4. Elimination of accumulative rounding error hy reorthogonalization process used to 
fit spherical harmonic functions of latitude and longitude to Fourier {time series) coefficients 
ai for UY2 monthly medians 

December 1957 113 stations. 





ORTHOGONAL 


COEFFiCiENTS 




DEVIATION OF 


ESTIMATES OF 


ACCUMULATIVE 












a^ IN PERCENT 


ROUNDING 


ERROR from: 


FIRST 




SECOND 




FIRST 


SECOND 




ORTHOGONALIZATION 


ORTHOGONALIZATION 




ORTHOGONALIZATION 


ORTHOGONALIZATION 












\-% 


AOk 


Abk 


k 


Ok 




\ 




bk 


Ok 


bk 





2.6578404E 


00 


2.6578404E 


00 


0. 


2.9041629E-06 


0. 


1 


1.5929625E 


00 


1.5929624E 


00 


-0.0000063 


2.6005150E-06 


0. 


2 


-1.4231348E 


00 


-1.4231393E 


00 


0.00032 


1.9663900E-05 


5.2353136E-08 


3 


-3.8013159E 


00 


-3.8013174E 


00 


0.000040 


-4.8294410E-06 


3.9199992E-08 


4 


-1.1098262E 


01 


-1.1098265E 


01 


0.000027 


7.0032999E-05 


0. 


5 


-1.4189972E 


01 


-1.4190128E 


01 


0.0011 


-1.1172427E-04 


8.4008608E-09 


6 


-1.0993559E 


01 


-1.0990924E 


01 


-0.024 


5.3348242E-03 


1.0846157E-08 


7 


-3.5136985E 


00 


-3.5109766E 


00 


-0.078 


-2.8899443E-02 


-2.5464985E-07 


8 


-1.7760643E 


00 


-1.8776263E 


00 


5.4 


1.1386564E 00 


6.2695740E-07 


9 


2.3313985E 


01 


2.3232524E 


01 


-0.35 


1.2002318E-02 


4.1049105E-08 


10 


1.7978380E 


02 


1.8443359E 


02 


2.5 


1.6151702E-02 


0. 


11 


1.0436464E- 


-01 


1.1301300E- 


-01 


7.7 


3.0752875E-02 


0. 


12 


-1.29641C1E- 


-01 


-1.2419979E- 


-01 


-4.4 


-3.1891573E-02 


7.4985841E-09 


13 


3.1378034E- 


-03 


3.4149244E- 


-03 


8.1 


3.2888361E-01 


7.6702861E-07 


14 


3.5222762E- 


-01 


3.3267845E- 


-01 


-5.9 


-3.2077405E-02 


1.1197871E-07 


15 


2.7970650E- 


-01 


2.4467795E- 


-01 


-14. 


-3.0703257E-02 


-2.0554128E-07 


16 


3.6830783E 


00 


3.6748188E 


00 


-0.22 


-5.6580728E-03 


0. 


17 


2.3149930E- 


-01 


1.7777535E- 


-01 


-30. 


-1.3896600E-01 


3.9814584E-07 


18 


-7.6929165E- 


-01 


-7.4671454E- 


-01 


-3.0 


-5.5304769E-02 


-1.4966725E-07 


19 


^4.8753538E 


00 


-4.7861179E 


00 


-1.9 


-3.4774151E-02 


3.7360954E-08 


20 


1.1794396E 


00 


1.0771025E 


00 


-9.5 


-1.6933961E-01 


-1.3834487E-07 


21 


7.3562072E 


00 


7,2621238E 


00 


-1.3 


-3.2170044E-03 


-3.2830421E-08 


22 


-1.5848168E 


00 


-1.3387020E 


00 


-18. 


-9.4834993E-02 


-7.7917363E-08 


23 


-1.2521099E 


01 


-1.2175401E 


01 


-2.8 


4.4063769E-03 


4.8954975E-08 


24 


-2.36B2257E 


01 


-2.3369081E 


01 


-1.3 


-3.2909461E-02 


-1.0202309E-08 


25 


8.9367033E- 


-02 


9.1821058E- 


-02 


2.7 


-6.0492957E-03 


-6.0856796E-08 


26 


5.9681830E- 


-01 


5.9769819E 


-01 


0.15 


1.3858429E-03 


0. 


27 


1.0897222E 


00 


1,1006036E 


00 


0.99 


1.1943727E-02 


0. 


28 


-2.7559001E 


-01 


-2.7433266E 


-01 


-0.46 


6.7750021E-03 


-2.7158927 -08 


29 


2.2791040E 


00 


2.2946516E 


00 


0.68 


6.1824901E-03 


2.5975466 -08 


30 


1.9115768E 


00 


1.9063635E 


00 


-0.27 


1.0157171E-04 


4.6899223 -08 
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B,ix)=A„{x) 



k-l 



where 



p=0 



bkn=- 



(A„B,) 



{Bp,B„) 



(50) 



(51) 



Then the least squares solution can be written in 
the form 

YK{x)=i:b,B,{x) (52) 

where 

I _ (y,B,) . . 

To obtain Fx(x) in the shnpler form (1), we compute 
first the coefficients al for expressing (52) in terms of 
the Aj^{x). This is clone by means of (13) and (14) 
with ajcy ajcp, and Djc replaced by 6^, bkp, and al, 
respectively. Then tlie desired coefficients D^ are 
obtained using (13) and (14) with Gjc replaced by al. 
The normalization of the solution is made by means 
of a similar modification of section 2.3. 

For illustration we continue with the same example 
given at the end of section 6.1. The coefficients 6a- 
and estimates of the relative error (Abjc/bjc) from one 
reorthogonalization are given in table 4. From 
these results it is seen that the largest value of the 
relative error is of the order of 10"^. Thus a great 
improvement lias been gained. Using tlie coeffi- 
cients bjc we have computed the relative deviation of 
the ajc in percent, from which it can be seen that 
an differs by as much as 30 percent. Moreover, 
the breakdown of the first orthogonalization at k=S 
is clearly shown. 

7. Summary of Analysis'^and" Numerical 
Maps 

We summarize here the steps used in the analysis 
described in the preceding sections (see table 5). 



The analysis begins with the actual observations as 
they are tabulated at the stations, each in its own 
zone time. The diurnal representation is obtained 
from Fourier analysis of the 24 hourly values (for 
each station), the corrections to LMT being pro- 
duced by an appropriate shift of the phase, det(\r- 
mined by the station locations. The optimum 
separation of the noise from the real physical vai'ia- 
tion of the data is made by truncating the liigh 
frequency harmonics (low-pass filtering). Thus for 
the characteristic /0F2, only 8 harmonics are needed. 
The geographic variation of each of the 17 Fourier 
coefficients 



aobiaib2a2 



(54) 



is represented by an orthonormal series, analogous 
to spherical harmonic analysis. Since the functions 
in these series must be orthogonal relative to the 
coordinates (X^, dt) of the stations, they cannot be 
known a priori; hence they nmst be constructed. 
Gram-Schmidt orthogonalization and reorthogonali- 
zation are used for this purpose. The smoothing 
(filtering) in the geographic dimensions is then per- 
formed by truncating the orthonormal series for 
each of the 17 Fourier coefficients. 

The end product of the analysis is a table of numer- 
ical coeffici(Mits D,„ dc^fining a function r(X,^,/) of 
the form 



T{\,d,t)=a,{\6)+J2[aj{\d) cos jt + bj{\,e) sinjf] (t>5) 
i=i 

whiTe each of the functions a^(X,^) and bj(\,d) — ^rep- 
resc^nting the geograpliic variations of the Foin^ier 
coefficients — is a series of the form 



k = 



(56) 



(see table 1). The index s denotes which Fourier 
coefficient is represented in the order shown in (54). 



Table 5. Summary diagram of the analysts 
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Table 6. Coefficients T>sk defining the function T (X, d, t) for monthly ynedian WF2 (Mc/s), December 1957 

1— main latitudinal variation. Mixed latitudinal variation: II— first order in longitude, III— second order in longitude. Notation: For each entry the number given by the first eight digits and sign is mul- 
tiplied by the power of 10 defined by the last two digits and sign. 

TIME VARIATION 



Harmonic 





1 


2 


3 


4 




^ 





1 


2 


3 


4 


5 


6 


7 


8 







1.0914267E 01 


7,7372089E-01 


7.283*942E-01 


-9.2623643E-01 


8.7422963E-01 


-1.7848102E-01 


-3.7885388E-01 


2,9685747E-01 


-2,5333249E-01 






-7.8953766E-01 


4.8827410E-01 


2,6530059E 00 


-1.5247263E 00 


8.0915360E-01 


-9.4284316E-01 


-1,7684497E 00 


-4.3241333E-01 


-3,7663186E-01 






2.9886163E 00 


1.6714562E 01 


2.1789511E 01 


6.5370734E-01 


-3.3546872E 00 


-2.3955230E-01 


-1.4676821E 00 


-1.4250625E 00 


-4,1630597E-01 






1.5785872E 01 


1.7292880E 00 


-1.8330226E 01 


3.4795284E 00 


-4,3719376E 00 


6.4229818E 00 


1.0586009E 01 


9.2466606E 00 


2,8930649E-02 






-7.3081383E 01 


-1.0531759E 02 


-1.0238298E 02 


2.1601855E 00 


8.6228765E CO 


-6.8914741E-01 


1.0163215E 01 


1.0430091E 01 


6.1317524E 00 


I 




-4.6525828E 01 


-4.1977129E 00 


7.5444276E 01 


2.0150497E 00 


1.6324815E 01 


-3.4641695E 01 


-2.6678432E 01 


-3.7404618E 01 


-5,6494755E-01 






1.9878676E 02 


2.5449798E 02 


2.4343033E 02 


-4.3352216E 00 


-5.4578'}38E 00 


-7.5223076E 00 


-2.2062303E 01 


-3.9342693E 01 


-2,0193537E 01 






4,9107151E 01 


8.4645829E 00 


-9,8828049E 01 


-2.6254765E 00 


-1,7752290E 01 


6,0494795E 01 


3.1563795E 01 


5.0754311E 01 


4.2224289E 00 






-2.1057577E 02 


-2.6631189E 02 


-2.6629260E 02 


1.1159760E 01 


-1.8912893E-01 


2.3693942E 01 


2.5115734E 01 


5.3091602E 01 


2.7875168E 01 






-1,7550868E 01 


-6.6770028E 00 


3.9049190E 01 


-1.4319514E 00 


5.0474002E 00 


-3.1529647E 01 


-1.37i8612E 01 


-2.2160749E 01 


-3.1592644E 00 




10 


7,6640091E 01 


1.0001173E 02 


1.0270474E 02 


-8.8129838E 00 


-5,9311831E-01 


-1.5101959E 01 


-1.1473124E 01 


-2.3054411E 01 


-1.3278131E 01 






-5.2964947E-01 


2.4364768E-03 


8.0027715E-01 


6.0721534E-01 


-1.2101191E-01 


-1.0403821E-01 


1.6257809E-01 


-1.0643796E-01 


6.2557100E-02 






2.6902534E-01 


9.4759300E-02 


-6.8859721E-01 


7,3459486E-01 


1,7502589E-01 


-5,0371087E-02 


3,1257521E-01 


-7.8402439E-02 


6.3306112E-02 






2.1349603E-02 


-3.9538481c 00 


3.4758507E-01 


3,0601999E 00 


-1,2871614E-01 


-4.8954588E-02 


9.1777242E-01 


-4.7842801E-01 


-1.8609319E-02 






3.1711371E 00 


1.7516087E 00 


-1.5318519E-01 


1.2964863E 00 


7,3990802E-01 


-8*4662998E-01 


2.2758654E-01 


-2.2977208E-02 


-1.0258091E-01 






1.6633779E 00 


-3.5908690E 00 


-1.1218609E 01 


-5.7828038E 00 


-5.0241930E-01 


7.0223299E-01 


-1.4589699E 00 


1.1074558E 00 


-5.0288115E-01 






-5.3166566E 00 


-5.6115813E 00 


-4,2479382E 00 


-7.8961277E 00 


7.7927373E-02 


2,7583131E 00 


-1.7414663E 00 


7.5675662E-01 


-8#2908363E-01 


II 




-1.91916^3E-01 


1.6253567E 01 


-4.2751713E 00 


-1.0306968E 01 


4.1552535E-01 


7,0501295E-01 


-4,5479047E 00 


8.5073202E-01 


-2.4932280E-01 




-2.?335084E 00 


-6.0404139E 00 


9,8574758E-02 


-6.3966695E 00 


-1.4409493E 00 


3.3823205E 00 


-1.9706686E 00 


1.9771550E-02 


-6.9674328E-03 






Z.1297312E-01 


8.3133165E 00 


3.0209412E 01 


1,2918815E 01 


2,0064190E 00 


-1.4150307E 00 


1,6666640E 00 


-1.5010523E 00 


8,8937262E-01 






4#1376244E 00 


2.0858187E 01 


2.1029316E 01 


1.8395100E 01 


-4.3110338E 00 


-1.2203444E 01 


2.3400389E 00 


-1.0840088E GO 


1.2606628E 00 








-1.346C412E 01 


5,3866998E 00 


7.1701701E 00 


-8,9934726E-02 


-1.3007687E 00 


4,4833956E 00 












3.7958689E 00 


2.9883151E-01 


6.3709797E 00 


9.0028155E-01 


-3#2107806E 00 


1.9325332E 00 












-6.2877868E 00 


-2.2921706E 01 


-8.3322061E 00 


-2.3156345E-01 


1.1989445E 00 


6,0279800E-01 












-1.6428649E 01 


-1.6010904E 01 


-1.1403061E 01 


4,9586704E 00 


1.1153906E 01 


-6.5103689E-01 












4.6818604E-01 


-8.7875202E-01 


-1.8927348E-02 


2.8565115E-01 


8,5442873E-02 


-3,1833852E-02 












-1.4529204E-01 


5.3212479E-01 


1.5694889E-01 


-U0888240E-01 


1.7714674E-03 


1.4413557E-02 






III 






-9.6026849E-02 


9.0155041E-01 


2.4230455E-01 


-5.1621379E-02 










28 




-2.7047523E-01 


-7.5395824E-01 


4.0790479E-01 


8,3877255E-02 












29 




-1.4390941E 00 


1.8835552E 00 
















30 




1.3630296E 00 


1.9063635E 00 
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Harmonic 


5 


6 


7 


8 




>^ 


9 10 


11 12 


13 14 


15 16 


I 




1 

2 
3 
4 
5 
6 
7 
8 
9 
10 


1.5531158E-01 2.3792567E-01 

3.7190015E-01 -4.9630663E-01 

1.0722732E-01 -6 •4481278E-01 

-2.7126487E-01 7,9842262E 00 

-2.0074097E 00 -9,5196921E-01 

3.0534070E-01 -3,1520131E 01 

6,2490299E 00 1 .6007400E-01 

-2.5598377E 00 4.3334295E 01 

-8.5661732E 00 4.7626390E 00 

2.1710776E 00 -1,9401572E 01 

4.0764945E 00 -3.6110235E 00 


-6.4560618E-02 6.1802934E-02 

3.7675107E-01 3 • 1683339E-01 

1.4638798E 00 -5 .7878543E-01 

-1.8491821E 00 1 . 1859947E-01 

-1.0122759E 01 6.0368895E 00 

3.9806334E 00 -2,7641763E 00 

2.9222164E 01 -1,9459016E 01 

-3.4801114E 00 3,3851973E 00 

-3,4945636E 01 2.3111368E 01 

9.0533099E-01 -1.0657850E 00 

1.4435670E 01 -9.1508738E 00 


-7.2905440E-02 -6. 1891910E-02 

-1,4461850E-01 2. 4323151E-01 

-3.7733971E-01 -3,0163709E-01 

-1.4023915E 00 -2,1925556E 00 

3,5854293E 00 4,5672736E 00 

7,4862415E 00 7.7719395E 00 

-8.4682503E 00 -1,1667597E 01 

-1.0514823E 01 -1.0894292E 01 

7,8529029E 00 1.1465121E 01 

4,6367167E 00 5.1017995E 00 

-2.4935791E 00 -4,0293859E 00 


7.4823447E-02 -1 •3109944E-03 

-1.0825114E-01 -1.0503997E-01 

-5.3879264E-01 -9,2140743E-03 

7.0639567E-01 -3. 1251254E-01 

1.4057245E 00 -1.6983739E 00 

-2.5793256E 00 1,9127105E 00 

-3.4862643E 00 7.5111334E 00 

3.9185429E 00 -2,6382915E 00 

5,0508122E 00 -1.0165159E 01 

-1.9346376E 00 1.1577847E 00 

-2,4991334E 00 4.3718649E 00 


II 


11 
12 
13 
14 
15 
16 


5,2536049E-03 4 • 544 1 1 54E-03 
-5.3775538E-02 -3 . 1418646E-02 

7.1804160E-02 5 ,29922 50E-02 

6.1635942E-02 -3 .0900826E-02 
-1.1925527E-02 -5 . 1757125E-02 

1.0552388E-01 4.9768387E-02 


2,7549960E-02 4 •8208075E-02 

6.6945602E-02 2 .2478048E-02 

2.2264430E-02 2 .715953 lE-02 

-8,3197396E-02 5 ,474276 lE-02 

-5.0750350E-02 -7. 1205906E-02 

-1.1215977E-01 -6#370^289E-02 


9.6875999E-04 -8. 5190815E-03 

-2.8802577E-02 1. 0904388E-02 

1.7879777E-03 3.3803478E-03 

3.2164980E-02 -1.8425275E-03 


-5.1576599E-03 -8 . 1442986E-04 
-2.6393688E-02 1 .2157885E-02 
-4,6019688E-03 -2 •2794990E-02 
4.8023952E-02 1 •0659245E-03 
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The runction T{\,d;t) tlierel'ore represents the con- 
tinuous time variations of an ionospheric charac- 
teristic on a worldwide basis and can be used to com- 
pute its value at any desired location and instant of 
time. Such a function is rc^ferred to as a numerical 
map. The coefficients defining the numerical map of 
foF2 median for December 1957 are given in table 6. 
Many useful applications are made from numerical 
maps. For example, a variety of worldwide contom* 
maps and other graphical representations of iono- 
spheric characteristics can be computed and plotted 
automatically by machine. As an illustration, 
figures 12 and 13 show the variation of /0F2 median 
for fixed instants of universal time (UT) and LMT, 
respectively. In addition, a great deal can be and 
has been learned about the regularities of the iono- 
sphere from numerical representations. The most 
important immediate application, however, is the 
prediction of long-term clianges in ionospheric con- 
ditions. Such predictions have already been pro- 
duced using the m(^thods described here. A further 
discussion of these results and otlier applications of 
numei'ical mapping will be given in later papers. 
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Figure 12. Maps of inonthly median foF2 (Mc/8) 

December 1957 



in U7iiversal time front the funelion T (X, 0, t). 
UT=12h. 
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Figure 13. Maps of monthly median foF2 (Mc/s) in local mean lime from the function T (X, d, t). 

December 1957 LMT=12h. 
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